home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
matrix.lha
/
Matrix
/
src
/
Matrix.ccP
< prev
next >
Wrap
Text File
|
1993-07-07
|
20KB
|
902 lines
pushdef(`index', `ifelse($#, 0, ``index'', ``index($*)'')')dnl
pushdef(`shift', `ifelse($#, 0, ``shift'', ``shift($*)'')')dnl
pushdef(`format', `ifelse($#, 0, ``format'', ``format($*)'')')dnl
pushdef(`index', `ifelse($#, 0, ``index'', ``index($*)'')')dnl
pushdef(`shift', `ifelse($#, 0, ``shift'', ``shift($*)'')')dnl
pushdef(`format', `ifelse($#, 0, ``format'', ``format($*)'')')dnl
#include <<T>.Matrix.h>
<T>Matrix <T>Array::operator - () const {
<T>* newX = new <T> [M*N];
<T>* t = newX;
<T>* u = X;
<T>* v = X;
do
while (u < v + N)
*t++ = -(*u++);
while ((u = v += L) < &X[M*L]);
return <T>Matrix(M, N, newX);
}
dnl 1$ = {*,/,+,-}
define(arithmetic,
`<T>Matrix operator $1 (const <T&> a, const <T>Array& b) {
<T>* newx = new <T> [b.m()*b.n()];
<T>* t = newx;
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v)
*t++ = a $1 *u++;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return <T>Matrix(b.m(), b.n(), newx);
}
<T>Matrix operator $1 (const <T>Array& a, const <T&> b) {
<T>* newx = new <T> [a.m()*a.n()];
<T>* t = newx;
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v)
*t++ = *u++ $1 b;
u += a.l() - a.n();
} while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
return <T>Matrix(a.m(), a.n(), newx);
}
<T>Matrix operator $1 (const <T>Array& a, const <T>Array& b) {
if (a.n() == b.n()) {
if (a.m() == b.m()) { // a.n() == b.n()
<T>* newx = new <T> [b.m()*b.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
<T>* v = b.x();
while ((t += b.n()) <= &newx[b.m()*b.n()]) {
while (s < t)
*s++ = *u++ $1 *v++;
u += a.l() - a.n();
v += b.l() - b.n();
};
return <T>Matrix(b.m(), b.n(), newx);
};
if (a.m() == 1) { // a.n() == b.n() && a.m() != b.m()
<T>* newx = new <T> [b.m()*b.n()];
<T>* s = newx;
<T>* t = newx;
<T>* v = b.x();
while ((t += b.n()) <= &newx[b.m()*b.n()]) {
<T>* u = a.x();
while (s < t)
*s++ = *u++ $1 *v++;
v += b.l() - b.n();
};
return <T>Matrix(b.m(), b.n(), newx);
};
if (b.m() == 1) { // a.n() == b.n() && a.m() != b.m() && a.m() != 1
<T>* newx = new <T> [a.m()*a.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
while ((t += a.n()) <= &newx[a.m()*a.n()]) {
<T>* v = b.x();
while (s < t)
*s++ = *u++ $1 *v++;
u += a.l() - a.n();
};
return <T>Matrix(a.m(), a.n(), newx);
};
};
if (a.n() == 1) { // b.n() != a.n()
if (a.m() == b.m()) { // b.n() != a.n() == 1
<T>* newx = new <T> [b.m()*b.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
<T>* v = b.x();
while ((t += b.n()) <= &newx[b.m()*b.n()]) {
while (s < t)
*s++ = *u $1 *v++;
u++;
v += b.l() - b.n();
};
return <T>Matrix(b.m(), b.n(), newx);
};
if (a.m() == 1) { // b.n() != a.n() == 1 && a.m() != b.m()
<T>* newx = new <T> [b.m()*b.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
<T>* v = b.x();
while ((t += b.n()) <= &newx[b.m()*b.n()]) {
while (s < t)
*s++ = *u $1 *v++;
v += b.l() - b.n();
};
return <T>Matrix(b.m(), b.n(), newx);
};
};
if (b.n() == 1) { // a.n() != b.n() && a.n() != 1
if (a.m() == b.m()) { // a.n() != b.n() == 1
<T>* newx = new <T> [a.m()*a.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
<T>* v = b.x();
while ((t += a.n()) <= &newx[a.m()*a.n()]) {
while (s < t)
*s++ = *u++ $1 *v;
u += a.l() - a.n();
v++;
};
return <T>Matrix(a.m(), a.n(), newx);
};
if (b.m() == 1) { // a.n() != b.n() == 1 && a.m() != b.m()
<T>* newx = new <T> [a.m()*a.n()];
<T>* s = newx;
<T>* t = newx;
<T>* u = a.x();
<T>* v = b.x();
while ((t += a.n()) <= &newx[a.m()*a.n()]) {
while (s < t)
*s++ = *u++ $1 *v;
u += a.l() - a.n();
};
return <T>Matrix(a.m(), a.n(), newx);
};
};
a.error("nonconformant <T>Array $1 operands.");
return <T>Matrix();
}
')dnl
arithmetic(*)
arithmetic(/)
<T>Matrix operator % (const <T>Array& a, const <T>Array& b) { // inner product
if (a.n() != b.n())
a.error("nonconformant <T>Array % operands.");
<T>* newx = new <T> [a.m()*b.m()];
<T>* s = newx;
<T>* t = a.x();
<T>* u;
<T>* v;
do {
v = b.x();
while (v < &b.x()[b.m()*b.l()]) {
u = t; *s = *u++ * *v++;
while (u < t + a.n())
*s += *u++ * *v++;
++s;
v += b.l() - b.n();
};
} while((t += a.l()) < &a.x()[a.m()*a.l()]);
return <T>Matrix(a.m(), b.m(), newx);
}
arithmetic(+)
arithmetic(-)
ifelse(basetype, complex, `
extern int Array_number; extern char* Array_format; extern char* Array_space;
', `
int Array_number = 4;
char* Array_format = "%g";
char* Array_space = " ";
char empty[1] = "";
char* format(char* s, int n, char* w) {
Array_number = n;
Array_format = s;
Array_space = w;
return empty;
}
')dnl
ostream& operator << (ostream& s, const <T>Array& b) {
<T>* t = b.x();
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
do {
t += (0 < Array_number)? Array_number: b.n();
t = (t < v)? t: v;
while (u < t)
if (!(ifelse(basetype, complex,
`s << "(" && s << form(Array_format, real(*u)) &&
s << ", " && s << form(Array_format, imag(*u)) &&
s << ")"',
`s << form(Array_format, *u)')
&& ((++u < t) ? s << Array_space : s << "\n")))
return s;
} while (t < v);
t = u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return s;
}
istream& operator >> (istream& s, const <T>Array& b) {
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v && s >> *u)
u++;
if (u < v)
return s;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return s;
}
<T>Matrix operator << (const <T>Array& a, int n) {
<T>* newx = new <T> [a.m()*a.n()];
<T>* t = a.x() + n;
<T>* u = newx;
<T>* v = newx;
while ((v += a.n()) <= &newx[a.m()*a.n()]) {
while (u < v - n)
*u++ = *t++;
while (u < v)
*u++ = 0.0;
t += a.l() - a.n() + n;
};
return <T>Matrix(a.m(), a.n(), newx);
}
<T>Matrix operator >> (const <T>Array& a, int n) {
<T>* newx = new <T> [a.m()*a.n()];
<T>* t = &a.x()[a.m()*a.l()];
<T>* u = &newx[a.m()*a.n()];
<T>* v = &newx[a.m()*a.n()];
while (newx <= (u -= a.n())) {
t -= a.l() - a.n() + n;
while (v > u + n)
*--v = *--t;
while (v > u)
*--v = *t;
};
return <T>Matrix(a.m(), a.n(), newx);
}
dnl 1$ = {==,< }
define(comparison,
`int operator $1 (const <T>Array& a, const <T&> b) {
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v && *u $1 b)
u++;
if (u < v)
return 0;
u += a.l() - a.n();
} while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
return 1;
}
int operator $1 (const <T>Array& a, const <T>Array& b) {
if (a.n() == b.n()) {
if (a.m() == b.m()) { // a.n() == b.n()
<T>* t = a.x();
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v && *t $1 *u) {
t++;
u++;
};
if (u < v)
return 0;
t += a.l() - a.n();
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return 1;
};
if (a.m() == 1) { // a.n() == b.n() && a.m() != b.m()
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
<T>* t = a.x();
while (u < v && *t $1 *u) {
t++;
u++;
};
if (u < v)
return 0;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return 1;
};
if (b.m() == 1) { // a.n() == b.n() && a.m() != b.m() && a.m() != 1
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
<T>* t = b.x();
while (u < v && *u $1 *t) {
t++;
u++;
};
if (u < v)
return 0;
t += b.l() - b.n();
u += a.l() - a.n();
} while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
return 1;
};
};
if (a.n() == 1) { // b.n() != a.n()
if (a.m() == b.m()) { // b.n() != a.n() == 1
<T>* t = a.x();
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v && *t $1 *u)
u++;
if (u < v)
return 0;
t++;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return 1;
};
if (a.m() == 1) { // b.n() != a.n() == 1 && a.m() != b.m()
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v && *a.x() $1 *u)
u++;
if (u < v)
return 0;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return 1;
};
};
if (b.n() == 1) { // a.n() != b.n() && a.n() != 1
if (a.m() == b.m()) { // a.n() != b.n() == 1
<T>* t = b.x();
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v && *u $1 *t)
u++;
if (u < v)
return 0;
t++;
u += a.l() - a.n();
} while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
return 1;
};
if (b.m() == 1) { // a.n() != b.n() == 1 && a.m() != b.m()
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v && *u $1 *b.x())
u++;
if (u < v)
return 0;
u += a.l() - a.n();
} while ((v += a.l()) <= &a.x()[a.m()*a.l()]);
return 1;
};
};
a.error("nonconformant <T>Array $1 operands.");
return 0;
}
')dnl
ifelse(basetype, complex, `',
`int operator < (const <T&> a, const <T>Array& b) {
<T>* u = b.x();
<T>* v = b.x() + b.n();
do {
while (u < v && a < *u)
u++;
if (u < v)
return 0;
u += b.l() - b.n();
} while ((v += b.l()) <= &b.x()[b.m()*b.l()]);
return 1;
}
comparison(`< ')
')dnl
comparison(==)
<T>Matrix operator & (const <T>Array& a, const <T>Array& b) { // outer product
if (a.m() != b.m())
a.error("nonconformant <T>Array & operands.");
<T>* newx = new <T> [a.n()*b.n()];
<T>* p = newx;
<T>* s = a.x();
<T>* t;
<T>* u;
<T>* v;
do {
t = b.x();
do {
u = s; v = t; *p = *u * *v;
while ((u += a.l()), (v += b.l()) < &b.x()[b.m()*b.l()])
*p += *u * *v;
p++;
} while (++t < &b.x()[b.n()]);
} while (++s < &a.x()[a.n()]);
return <T>Matrix(a.n(), b.n(), newx);
}
<T>Matrix operator ^ (const <T>Array& a, const <T>Array& b) { // stack
if (a.n() != b.n())
a.error("nonconformant <T>Array ^ operands.");
<T>* newx = new <T> [(a.m()+b.m())*a.n()];
<T>* t = newx;
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v)
*t++ = *u++;
u += a.l() - a.n();
} while((v += a.l()) <= &a.x()[a.m()*a.l()]);
u = b.x();
v = b.x() + b.n();
do {
while (u < v)
*t++ = *u++;
u += b.l() - b.n();
} while((v += b.l()) <= &b.x()[b.m()*b.l()]);
return <T>Matrix(a.m()+b.m(), a.n(), newx);
}
<T>Matrix operator | (const <T>Array& a, const <T>Array& b) { // adjoin
if (a.m() != b.m())
a.error("nonconformant <T>Array | operands.");
<T>* newx = new <T> [a.m()*(a.n()+b.n())];
<T>* t = newx;
<T>* u = a.x();
<T>* v = a.x() + a.n();
do {
while (u < v)
*t++ = *u++;
t += b.n();
u += a.l() - a.n();
} while((v += a.l()) <= &a.x()[a.m()*a.l()]);
t = newx + a.n();
u = b.x();
v = b.x() + b.n();
do {
while (u < v)
*t++ = *u++;
u += b.l() - b.n();
t += a.n();
} while((v += b.l()) <= &b.x()[b.m()*b.l()]);
return <T>Matrix(a.m(), a.n()+b.n(), newx);
}
dnl 1$ = { ,*,/,+,-}
define(assignment,
`<T>Array& <T>Array::operator $1= (const <T&> b) {
<T>* s = X;
<T>* t = X + N;
do {
while (s < t)
*s++ $1= b;
s += L - N;
} while ((t += L) <= &X[M*L]);
return *this;
}
<T>Array& <T>Array::operator $1= (const <T>Array& a) {
if (N == a.N) {
if (M == a.M) { // N == a.N
<T>* t = a.X;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*u++ $1= *t++;
t += a.L - a.N;
u += L - N;
} while ((v += L) <= &X[M*L]);
return *this;
};
if (a.M == 1) { // N == a.N && M != a.M
<T>* t;
<T>* u = X;
<T>* v = X + N;
do {
t = a.X;
while (u < v)
*u++ $1= *t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return *this;
};
};
if (a.N == 1) { // N != a.N
if (M == a.M) { // N != a.N == 1
<T>* t = a.X;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*u++ $1= *t;
t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return *this;
};
if (a.M == 1) { // N != a.N == 1 && M != a.M
<T>* t = a.X;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*u++ $1= *t;
u += L - N;
} while ((v += L) <= &X[M*L]);
return *this;
};
};
error("nonconformant <T>Array $1= operands.");
return *this;
}
')dnl
assignment(` ')
assignment(*)
assignment(/)
assignment(+)
assignment(-)
<T>Array& <T>Array::operator >>= (int n) {
<T>* t = &X[M*L];
<T>* u = &X[M*L];
<T>* v = &X[M*L];
while (X <= (u -= L)) {
t -= L - N + n;
v -= L - N;
while (t > u)
*--v = *--t;
while (v > u)
*--v = *t;
};
return *this;
}
<T>Array& <T>Array::operator <<= (int n) {
<T>* t = X + n;
<T>* u = X;
<T>* v = X + N;
do {
while (t < v)
*u++ = *t++;
while (u < v)
*u++ = 0.0;
t += L - N + n;
u += L - N;
} while ((v += L) <= &X[M*L]);
return *this;
}
// Error Handling
void default_<T>Array_error_handler(const char* msg) {
cerr << "Fatal <T>Array error. " << msg << "\n";
exit(1);
return;
}
one_arg_error_handler_t
<T>Array_error_handler = default_<T>Array_error_handler;
one_arg_error_handler_t set_<T>Array_error_handler(one_arg_error_handler_t f) {
one_arg_error_handler_t old = <T>Array_error_handler;
<T>Array_error_handler = f;
return old;
}
void <T>Array::error(const char* msg) const {
(*<T>Array_error_handler)(msg);
}
ifelse(basetype, complex, `', `
<T>Matrix <T>Array::i(double epsilon) const { // inverse
<T>* a = new <T> [N*M];
<T>* uu = new <T> [M*N];
<T>* vv = new <T> [N*N];
<T>* w = new <T> [N];
<T>** u = new <T>* [M];
<T>** v = new <T>* [N];
void svdcmp(<T>**, int, int, <T>*, <T>**);
for (int i = 0; i < M; i++)
u[i] = &(uu[N*i]);
for (int j = 0; j < N; j++)
v[j] = &(vv[N*j]);
for (i = 0; i < M; i++)
for (j = 0; j < N; j++)
u[i][j] = X[i*L+j];
svdcmp(u, M, N, w, v); // Singular value decomposition.
<T> wmax = 0.0; // Maximum singular value.
for (j = 0; j < N; j++)
if (w[j] > wmax)
wmax = w[j];
<T> wmin = wmax*epsilon;
for (int k = 0; k < N; k++)
if (w[k] < wmin)
w[k] = 0.0;
else
w[k] = 1.0/w[k];
for (i = 0; i < N; i++)
for (j = 0; j < M; j++) {
a[M*i+j] = 0.0;
for (k = 0; k < N; k++)
a[M*i+j] += v[i][k]*w[k]*u[j][k];
};
delete [] w;
delete [] u;
delete [] v;
delete [] uu;
delete [] vv;
return <T>Matrix(N, M, a);
}
')
<T>Matrix <T>Array::t() const { // transpose
<T>* newX = new <T> [N*M];
<T>* t = newX;
<T>* u;
<T>* v = X;
do {
u = v;
do
*t++ = *u;
while ((u += L) < &X[M*L]);
} while (++v < &X[N]);
return <T>Matrix(N, M, newX);
}
<T>Matrix <T>Array::sum() const {
<T>* newX = new <T> [M];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
*t = *u;
while (++u < v)
*t += *u;
t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(1, M, newX);
}
<T>Matrix <T>Array::sumsq() const {
<T>* newX = new <T> [M];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
*t = *u * *u;
while (++u < v)
*t += *u * *u;
t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(1, M, newX);
}
ifelse(basetype, complex, `
#ifdef __ATT_complex__
<T>Matrix <T>Array::map(const <T> (*f)(const <T>&)) const {
<T>* newX = new <T> [M*N];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*t++ = f(*u++);
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(M, N, newX);
}
#else
<T>Matrix <T>Array::map(const <T> (*f)( <T> )) const {
<T>* newX = new <T> [M*N];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*t++ = f(*u++);
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(M, N, newX);
}
#endif
doubleMatrix <T>Array::map(const double (*f)(const <T>&)) const {
double* newX = new double [M*N];
double* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*t++ = f(*u++);
u += L - N;
} while ((v += L) <= &X[M*L]);
return doubleMatrix(M, N, newX);
}
#ifndef __ATT_complex__
doubleMatrix <T>Array::map(const double (*f)( <T> )) const {
double* newX = new double [M*N];
double* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*t++ = f(*u++);
u += L - N;
} while ((v += L) <= &X[M*L]);
return doubleMatrix(M, N, newX);
}
#endif
',
`<T>Matrix <T>Array::map(const <T> (*f)( <T> )) const {
<T>* newX = new <T> [M*N];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
while (u < v)
*t++ = f(*u++);
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(M, N, newX);
}
<T>Matrix <T>Array::min() const {
<T>* newX = new <T> [M];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
*t = *u;
while (++u < v)
if (*t > *u)
*t = *u;
t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(1, M, newX);
}
<T>Matrix <T>Array::max() const {
<T>* newX = new <T> [M];
<T>* t = newX;
<T>* u = X;
<T>* v = X + N;
do {
*t = *u;
while (++u < v)
if (*t < *u)
*t = *u;
t++;
u += L - N;
} while ((v += L) <= &X[M*L]);
return <T>Matrix(1, M, newX);
}
int <T>Array::min_index() const {
<T>* t = X;
<T>* u = X;
<T>* v = X + N;
do {
do
if (*t > *u)
t = u;
while (++u < v);
u += L - N;
} while ((v += L) <= &X[M*L]);
return (t - X);
}
int <T>Array::max_index() const {
<T>* t = X;
<T>* u = X;
<T>* v = X + N;
do {
do
if (*t < *u)
t = u;
while (++u < v);
u += L - N;
} while ((v += L) <= &X[M*L]);
return (t - X);
}
')dnl
ifelse(basetype, complex, `
<T>Matrix fft(const <T>Array& a) {
extern void four1(double*, int, int);
<T>Matrix result(a);
for (int i = 0; i < result.m(); i++)
four1((double*) result[i] - 1, result.n(), 1);
return result;
}
<T>Matrix fft(const doubleArray& a) {
extern void realft(double*, int, int);
<T>Matrix data(a.m(), a.n());
doubleArray(2*a.n(), a.m(), a.n(), (double*) data.x()) = a;
for (int i = 0; i < a.m(); i++) {
realft((double*) data[i] - 1, a.n()/2, 1);
for (int j = 1; j < a.n()/2; j++)
data[i][a.n()-j] = conj(data[i][j]);
data[i][a.n()/2] = imag(data[i][0]);
data[i][0] = real(data[i][0]);
};
return data;
}
', `
<T>Matrix ffct(const <T>Array& a) {
extern void cosft(double*, int, int);
<T>Matrix y(a);
for (int i = 0; i < a.m(); i++)
cosft(y[i] - 1, a.n(), 1);
return y;
}
')
// Constructors
<T>Matrix::<T>Matrix(const <T>Matrix& a) :
<T>Array(a.M, a.N, new <T> [a.M*a.N]) {
<T>* t = a.X;
<T>* u = X;
<T>* v = X;
while ((v += N) <= &X[M*N]) {
while (u < v)
*u++ = *t++;
t += a.L - a.N;
};
}
<T>Matrix::<T>Matrix(const <T>Array& a) :
<T>Array(a.M, a.N, new <T> [a.M*a.N]) {
<T>* t = a.X;
<T>* u = X;
<T>* v = X;
while ((v += N) <= &X[M*N]) {
while (u < v)
*u++ = *t++;
t += a.L - a.N;
};
}
ifelse(basetype, complex, `
<T>Matrix::<T>Matrix(const doubleArray& r, const double i) :
<T>Array(r.M, r.N, new <T> [r.M*r.N]) {
<T>* u = X;
<T>* v = X;
double* t = r.X;
while ((v += N) <= &X[M*N]) {
while (u < v)
*u++ = complex(*t++, i);
t += r.L - r.N;
};
}
<T>Matrix::<T>Matrix(const doubleArray& r, const doubleArray& i) :
<T>Array(r.M, r.N, new <T> [r.M*r.N]) {
if(i.M == M && i.N == N) {
<T>* u = X;
<T>* v = X;
double* t = r.X;
double* k = i.X;
while ((v += N) <= &X[M*N]) {
while (u < v)
*u++ = complex(*t++, *k++);
t += r.L - r.N;
k += i.L - i.N;
};
}
else
error("nonconformant <T>Array operands.");
return;
}
')dnl
<T>Matrix::~<T>Matrix()
{ delete [] X; }
popdef(`index')dnl
popdef(`shift')dnl
popdef(`format')dnl